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We study superconductivity in an ultracold Bose-Fermi mixture loaded into a square optical lat- 
tice subjected to a staggered flux. While the bosons form a superfluid at very low temperature 
and weak interaction, the interacting fermions experience an additional long-ranged attractive in- 
teraction mediated by phonons in the bosonic superfluid. This leads us to consider a generalized 
Hubbard model with on-site and nearest-neighbor attractive interactions, which give rise to two 
competing superconducting channels. We use the Bardeen-Cooper-Schrieffer theory to determine 
the regimes where distinct superconducting ground states are stabilized, and flnd that the non-local 
pairing channel favors a superconducting ground state which breaks both the gauge and the lattice 
symmetries, thus realizing unconventional superconductivity. Furthermore, the particular structure 
of the single-particle spectrum leads to unexpected consequences, for example, a dome-shaped su- 
perconducting region in the temperature versus flling fraction phase diagram, with a normal phase 
that comprises much richer physics than a Fermi-liquid. Notably, the relevant temperature regime 
and coupling strength is readily accessible in state of the art experiments with ultracold trapped 
atoms. 



I. INTRODUCTION 

Superconductivity in low dimensional systems is more 
than ever an active field of research. Despite decades 
of research, one of its holy grail remains to be a thor- 
ough understanding of high-temperature superconduc- 
tivity (high-Tc) in the cupratcs [1, 2]. Even though the 
Bardeen-Cooper-Schrieffer (BCS) theory has been proven 
to be very successful in addressing ordinary supercon- 
ducting phenomena, its straightforward application to 
understand the strongly correlated regime (e.g. in high- 
Tc superconductors) remains a challenge. 

In this article we argue that, in combination with un- 
conventional single-particle spectra, BCS theory gives 
rise to phenomena reminiscent of the physics known 
to occur in strongly correlated systems. We therefore 
study a iiltracold atom system with a Dirac-likc spectrum 
with linear rather than the ordinary quadratic disper- 
sion, arising via a time-reversal symmetry braking term 
in the Hamiltonian. At first sight, the relevance of Dirac 
fermions appears to be limited to a relativistic context. 
However, such fermions do emerge in condensed matter 
systems under certain circumstances. Examples include 
high-Tc cuprate superconductors, which exhibit a d^^-y^ 
symmetry in the superconducting order parameter, such 
that fermionic excitations along the nodal lines on the 
Fermi surface are Dirac-like [3] , and graphene, where the 
hexagonal crystal lattice gives rise to Dirac-like excita- 
tions [4, 5]. 

Indeed, the recent breakthrough in fabricating sheets 
of graphene has provided us with a solid-state model of 
two-dimensional Dirac fermions [6]. The observation of 
the half-integer quantum Hall effect [7] and Klein tun- 
neling [8, 9] in graphene, for example, are hallmarks 
of two-dimensional relativistic physics taking place in a 



condensed matter system. Interesting theoretical works 
have explored the importance of interaction effects in 
graphene, where phenomena such as room-temperature 
superfluidity [10], an anomalously low shear viscosity in 
the vicinity of quantum criticality [11] and novel super- 
conductivity [12, 13] have been predicted. On the exper- 
imental front, the anticipated fractional quantum Hall 
effect has only been observed very recently in an exfoli- 
ated graphene sample [14, 15]. Nevertheless, it still re- 
quires ingenuity to prepare a clean, highly-controllable 
graphenc-based system, in order to fully explore its rich 
physics. 

In ultracold atomic systems, on the other hand, the 
remarkable progress made in the last decade has al- 
lowed for experimental demonstrations of prototypical 
many-body phenomena, as the superfluid-Mott insulator 
transition in the Bose-Hubbard model [16-19] and the 
crossover between the Bardeen-Cooper-Schrieffer (BCS)- 
Bose-Einstein condensation (BEC) regimes [20-25]. One 
of the major tasks in the field is the use of optical lattices 
for systematic emulations of the Hubbard model, a model 
which is believed to capture the physics of high-Tc su- 
perconductors. A particular interesting option in optical 
lattices is the study of degenerate Bose-Fermi mixtures. 
The relative ease in tuning the parameters such as the in- 
terspecies interaction strength, and the density and mass 
ratios, offers a wide window for studies confronting the- 
ory with precision experiments. Many interesting phe- 
nomena have been investigated, which include the pre- 
diction of exotic quantum matters, such as supersolids 
[26-28], composite fermions [29], charge density waves 
and polaron-like quasi-particles [30]. Also more com- 
mon effects arise, such as the enhancement of different 
types of fermionic superfluidity mediated by the phonon 
background provided by the bosons [31, 32], which would 
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FIG. 1: (color online) (a) Schematic of the two-dimensional 
square optical lattice with a spacing of half the laser wave- 
length A/2. Under the presence of a flux passing through 
each plaquette that alternates in sign across the neighboring 
plaquette - staggered flux - it is then convenient to split the 
lattice into two inequivalent sublattices A and B. (b) The 
Bravais lattice is made of the ^-sublattice with a lattice con- 
stant d = \/2A/2. The unit cell consists of a basis, which is 
defined on the Z3-sublattice. The four vectors e; , with I — 1—4, 
connect each ^-site to its four nearest-neighboring B-sites. 



otherwise be difficult to control in other systems. On the 
experimental side, the ability to control the experimental 
parameters in Bose-Fermi mixtures in optical lattices is 
steadily progressing [33-35]. 

In this work, we study the system of Bose-Fermi mix- 
tures in a two-dimensional optical square lattice that pro- 
vides an inherent staggered flux. In a recent letter we 
have shown the emergence of Dirac fermions and the sta- 
bilization of various superconducting states for this sys- 
tem [36] . Here, we present the details of the calculations 
and extend the work to include a discussion on the com- 
peting spin-liquid instability. In addition, we study the 
re-entrant behavior due to the specific pairing states and 
the evolution of the Fermi surface as a function of the 
staggered flux. 

The paper is organized as follows: in Sec. II, we use 
the Bogoliubov theory to study the mediation effects of 
the condensed bosons on the fermions in the presence of 
the staggered flux, thus providing a method to realize an 
extended Hubbard interaction. In Sec. Ill, we develop 
a mean-field theory to study superconductivity with the 
model. We first discuss possible pairing states in the 
separate interaction channels and their associated sym- 
metries. We then construct the free energy and obtain 
the gap equations. We find a re-entrant behavior, which 
is unique to the local pairing channel. As for the non- 
local pairing channel, we compare the transition temper- 
atures of the various possible pairing states. We then 
discuss the evolution of the Fermi surface as the stag- 
gered flux is tuned away from the case of a 7r-flux. In 
Sec. IV, experimental signatures are discussed for the 
detection of the various superconducting order parame- 
ters and the evolution of the Fermi surfaces. Finally, we 
close the paper with a discussion and conclusions in Sec. 
V. 



II. BOSE-FERMI MIXTURES 

We start from the microscopic single-band Bose-Fermi 
Hubbard Hamiltonian subjected to a staggered flux 
[37, 38], 
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Here, Jb and J are the hopping amplitudes for bosons 
and fermions, respectively, and Ubb, Uff and Ubf are 
the boson-boson, fermion-fermion and boson- fermion on- 
site interactions, respectively. The presence of the stag- 
gered magnetic field leads to the appearance of two in- 
equivalent ^- and S-sublattices, see Fig. 1(a). The oper- 
ators Or and br+ei are the bosonic annihilation operators 
acting on site r and r + ei of the A- and i3-sublattices, 
respectively, and a^^a and 6r+e,,(T are the corresponding 
fermionic annihilation operators with spin a. Finally, nf 
and the boson and fermion number operators on 

site r, respectively. 

We write the grand canonical partition function Z in 
the functional-integral representation 

Z = y" VaVa*Va^Valeyi^i^-^{SB + Sp + Si)^ (3) 

where the bosonic action Sb is given by 

M r 

Sb ^ dr Y a*{T){hdr - ^J.)ai{T) + Ho^B + Hbb 
the fermionic action 5"^^ is given by 
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and the boson-fermion interaction term Sj is given by 
Si — Hbf- 
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Here, /i is the chemical potential for bosons, and the 
one for fcrmions with spin a. The inverse thermal energy 
at temperature T is given by /3 = l/ksT. 

We focus on the weakly interacting regime of the 
atomic bosons where they form a Bose-Einstein conden- 
sate (BEC) in the lattice. Such regime can be described 
accurately within the Bogoliubov theory. Before we ap- 
ply the theory, we need to identify the condensation mode 
for the bosons, since the presence of the staggered flux 
can give rise to distinct BECs. 

We first note that the single-particle term iJo.s can 
be diagonalized in momentum space with the canonical 
transformation 
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is the lattice dispersion, and the new operators 
c^iciT) , /3i^{t) correspond to the upper and lower energy 
band states. The bosonic operator defined on the sublat- 
tice A can then be written as 
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where the lattice momentum k is defined in the first 
Brillouin zone {IBZ) and the canonical transformation 
is used. Note that the total number of lattice sites is 
A'' = 2Na- To simplify the following calculations, we 



map the upper band operator a\^{T) in the first Brillouin 
zone to the second Brillouin zone {2BZ) which is then 
denoted by /3k (''■)• Since the transformation coefficient 
changes an overall sign in the second Brillouin zone, we 
have 
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Similarly, for the bosonic operators on the B sublattice, 
we have 



^ik-(ri+ei)^ j^y-j 

kelBZ©2BZ 

We may now identify the condensation mode kg as the 
single-particle state with the lowest energy, and per- 
form the c-number substitution for the condensate field 
/3ko {t) ^ as follows 



^^giko.(r,+ei) ^ 
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no = Nq/N defines the condensate density. The prime in 
the momentum summation means that the k = kg term 
is omitted. 

By expanding also the fields in Matsubara frequencies 



'^k(Wm), 



(8) 



with ujm = 2Trm/hP, we make the Bogoliubov approxima- 
tion in the action Sb + Sj where the bosonic fiuctuation 
field /3k (t) is kept to the quadratic order to obtain (see 
Appendix) 
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where 

with the source term containing tlie fermionic fields given by 

I 



and E]^ = =F|ek| for k e 1BZ(2BZ), and 

M(ki,k2,k3,k4) EE 1 + .9ki.9k2 5k3 5k4- Several sym- 
metry properties have been employed: = g-k and 
ko = {±ko}, in the sense that for r = adi + one 
gets exp[iko • r] = exp[i7r(±a ± 7)]. We have chosen the 
bosonic chemical potential to obey the Hugenholtz-Pines 
form, which depends also on the fermion mean-field den- 
sity ha, 

11 = Eo + UBBno + UBF^h„, (12) 

cr 

SO that the bosonic fluctuation field remains massless (the 



Goldstone mode). The action Eq. (9) is now at most of 

quadratic order in the bosonic fluctuation fleld (j) and we 
can thus integrate it out analytically [39, 40]. 

Starting from the original action (3), by performing 
the Bogoliubov approximation and integrating out the 
bosons, we arrive at the effective action for the fermions, 

k6lBZ®2BZ m 
= Sp + Sind (13) 

where the induced interaction is given by 



k,m 

Here, short hand notations are used for the form factor M = M(k, — k, ko, kg), i^kko = (-^k — -^'ko + ^ss^o) aiid 
W^k = - |MpC/M/4. 

Next, we consider the frequency independent component of the induced interaction, a well-studied regime for Bose- 
Fermi mixtures (e.g., the widely used rubidium-potassium system [41]) in the lattice. Due to differences in the laser 
detuning as experienced by the different atomic species, the hopping amplitude can be realized for Jb ^ J so that 
the consideration of the static limit is justifled [31]. By setting huirn = and converting the momentum summation 
into a momentum integral (1/iV) J2\i ~^ d'^ Iibz®2BZ c^^k/(27r)^ we get 

Sint = ^ E ^^AA{TC-Y')[nr,anr',a' +nr+ei,anT'+e^,a']+ E ^^AB{r - r' - ex)nr,anr'+ei,,7' , (15) 

where the induced potentials are given by 

/d^k 2 f 
'{2^W^Y^'^° cos[ko • (r - r')] cos[k • (r - r')] 

-?7BBnoCOs[k ■ (r - r')] cos[ko • (r + r')] cos'lOn - ^k]| = (^k4, (16) 

and 
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(27r)2 Wk V*'"''" ^^^^ • (r - r' - ei)] cos[ko • (r - r' - ei) -9i, + Oq] 
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4C/|^^, /, Jb |r-r'-ei 



These are the nonlocal phonon-mediated fermion density- 



density interaction terms, which take the form of a 
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Yukawa potential (screened Coulomb) in two dimensions. 
Mediated by the phonons, a non-local attractive inter- 
action is generated between fermions of all spin states, 
which fall off on the scale of the healing length ^. Exper- 
iments in a bosonic 2D lattice of rubidium atoms by Spiel- 
man et al. [19] show that typical values of ^ are on the 
order of d/^/2. For Js/UBBno = 1, the induced poten- 
tials are ^423(61) ~ -(4C/|^/i7Bi3) 0.17 while VA^(di) 
is one order of magnitude smaller. 

Thus the effective interactions are the renormal- 
ized on-site interaction gi = —UpF — ^4^(0), and 
the nearest neighbor interaction 52 = ^y^B(ei) — 
{AU'^p/Ubb)^-^'^ ■ Note that the boson-boson interac- 
tion Ubb is taken to be positive for the stability of 
the Bose gas and thus, the induced interaction is al- 
ways attractive. We write the effective interaction for 
the fermionic Hamiltonian as 
_5i 
2 



Hi, 



^ ^ ^ ^ ^r,cr^r, — (T 
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Notice that since U pp can be tuned by changing the scat- 
tering length ttg via Feshbach resonances, the relative 
strength g\/g2 can be tuned in a straightforward man- 
ner. Furthermore, the dependence of the induced poten- 
tials on the parameter {Jb /Ubbit-q) can also be taken as 
an extra tuning parameter; both techniques thus allow 
for an independent control of both gi and g2- 



III. NOVEL SUPERCONDUCTIVITY 

In the last section, we studied the attractive mediation 
effect of the bosonic superfluid on the fermions, which 
is shown to extend to the nearest-neighbor sites. The 
effective Hamiltonian describing the interacting fermions 

H = Ho,F + (19) 

provides the basis for the study of superconducting in- 
stabilities due to the competing on-site and the nearest- 
neighbor attractive interactions. For clarity of presenta- 
tion, we first identify the various order parameters that 
are favored by the interaction terms and carry out the 
analysis for the different superconducting channels inde- 
pendently. Then, we construct the full phase diagram in 
the presence of both couplings. 



A. Pairing Hamiltonian 

For the on-site attractive interaction gi > 0, we con- 
sider an on-site spin-singlet pairing with the following 
superconducting order parameter 

ieA ieB 



Al 



N 



By performing a mean-field decoupling in Hi with re- 
spect to this superconducting order parameter and keep- 
ing fluctuations up to first order, we arrive at the follow- 
ing mean-field Hamiltonian 
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MF 



51 
-FH.c. 
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For the nearest-neighbor attractive interaction 52 > 0, 
we consider the following decomposition: 



<i,j> 



(21) 



The first term may give rise to a non-local pairing 
correlation involving two nearest-neighbor sites with a 
spin-singlet structure. It is the resonating-valence-bond 
(RVB) state that was proposed by P. W. Anderson [42] . 
We take the non-local pairing correlation to be of the 
form 

_2g2 
N 



A2(e() 



which results in a superconducting order parameter 



ike. 
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Note that this superconducting order parameter is deter- 
mined by four independent RVB components A2(ei). 
In the second and third terms of Eq. (21), a non-trivial 

expectation value amounts to a particle-hole correlation 

= (ar,t^i+e,,t + "••,4.^'i+ei,4.)- Corresponding order 
parameter is given by 



r = g2 
~ 2N 



ieA 1=1 



By performing a mean-field decoupling in the Hamilto- 
nian H2 with respect to the two orders, we obtain 
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where = = 4 cos(A:2;(i/2) cos(fcj^(i/2). It is 

worth noting that the contribution from the second-order 
particle-hole correlation is negative, as opposed to the 
contribution from the superconducting order, which can 
be traced back to an additional overall minus sign when 
performing the mean-field decoupling for the particle- 
hole terms in Eq. (21). 



B. Symmetry of the Order Parameters 

Before we discuss the self-consistent procedure to de- 
termine the mean-field state, it is important to under- 
stand the symmetry of the superconducting order param- 
eter in a basis where the single-particle term is diagonal. 
By transforming to that basis, the mean-field interaction 
terms become 



Hi, 



MF 



cos(6'k) (/3k,i/5-k,t + ak.itt-k,- 



H.c. 
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2,MF 



-isin(6lk)(/3k,ia-k,t + ak,;/3-k,t) 
^ A^ k e''" (/3k,4/3-k,t - ak,4a-k,t) + H.c. 



Here, ^k.u and ak,cr are the lower and upper band op- 
erators, which obey anti-commutation relations, 6^ = 
ek/|ek|, and the property Ck = e_k or 9k = 0-\s. follows 
from Eq. (5). In i?2,MF, we have omitted the particle- 
hole correlation contribution since we are only interested 
in the symmetry of the superconducting order parameter 
here. 

In the pairing potential Hi^mf, we first note that the 
superconducting order parameter Ai induces an intra- 

band pairing (/^Ikf^k^ + '^-kt^^ki-' ^® ^^^^ ^ 
inter-band pairing (/?lkt"ki '^-kt^^i) band 
representation. With the orbital part Aje'^'' being an 
even function in the momentum variable, the intra-band 
pairing naturally gives rise to a spin singlet structure. 
For the inter-band pairing, since it is symmetric in the 
band index, it also results in a spin singlet structure. 
Thus, the on-site pairing is described by a total order 
parameter that has a spin singlet structure. 

On the other hand, the RVB order A2{ei) induces only 
intra-band pairing (/3k,i/3-k.t " ak,iQ!-k,t) in H2,mf- 
The spin structure of the pairing is then determined by 
the parity of the superconducting order A2,k, which de- 
pends on the choice of the four-component RVB order. 
For an even parity A2,k — A2,_k, it results in a spin- 
singlet, whereas for an odd parity A2,k = — A2,-k. it 
results in a spin triplet. It is important to note that 
since the external staggered flux breaks the A-B sublat- 
tice symmetry explicitly, the order parameter A2,k need 
not have a definite parity. It can generally have a mixed 
parity with a coherent mixture of states with even and 
odd parities, and the resulting spin structure is also a 
coherent mixture of spin singlet and spin triplet states. 



C. Mean-Field Free Energy 

The mean-field procedure that we carried out until now 
involves the introduction of five superconducting order 
parameters, i.e. Ai and A2(ei), for / = 1,2,3,4, and 
a particle-hole correlation F. For self-consistency, we 
need to minimize also the resulting free energy. To eval- 
uate the free energy, we first rewrite the grand-canonical 
mean-field Hamiltonian in the following form 



Hmf = Eo+ J2 *kMk, 

kelBZ 
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and ©k is given by the following matrix 
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The fermionic chemical potential with equal spin pop- 
ulation is At) and the full lattice dispersion is given 
by Eq. (5). Since the mean- field Hamiltonian (23) 
is quadratic, we perform a canonical transformation 
(a Bogoliubov-Valatin transformation) to diagonalize it. 
The system is then described by two branches of non- 
interacting fermionic quasi-particles and quasi-holes with 
spectra ±£'y^k, with v = 1,2. In the new basis, the en- 
tropy of the system can be computed simply from that 
for the ideal Fermi gas 

S = -fcB^|/(£;,,k)ln/(£,,k) 
k 

+ [1 - /(^,.,k)] ln[l - ./(£.,k)] 

where f{E^^k) = l/[exp(^i?y^k) + 1] is the Fermi-Dirac 
distribution. The free energy is then given by 

F[Ai,A2(eO] = E-TS 



(24) 



' i/=l k 

+ ln(l + e^^-^) 



We see that the jFp-term renders the free energy un- 
bounded below, which is an artifact of the mean-field 
decoupling for a particle-hole correlation. Indeed, the 
mean-field procedure is not suited to treat such a corre- 
lation. We will therefore ignore its contribution and set 
it to zero F = for the rest of the work. In the following 
sections, we minimize the above expression with respect 
to the variational parameters (Ai,A2(e;)) and identify 
the regime of parameters where superconductivity can 
occur. 
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FIG. 2: (a) Solutions to the gap equation for local supercon- 
ductivity at the phase transition for different chemical poten- 
tials: /i/J = (solid), ^/J = O.f (dash), ^/J = f (dotted 
dash), /i/J = 2 (dotted), (b) To show the reentrant behavior, 
we plot the evolution of the free energy at different tempera- 
tures fcsT/J = 0.8 (dashed), 0.56 (solid), 0.1 (dotted), for a 
fixed coupling gi/J = 6.1 and chemical potential /i/J = 0. 



D. Local Superconductivity 

In this section we consider only the local pairing chan- 
nel characterized by the order parameter Ai. Setting 
^2,k = 0, the quasi-particle spectra can be readily ob- 
tained, 



E- 



l,2,k 



By substituting these into the free energy function 
Eq. (24) and extremizing the latter with respect to Ai, 



we obtain the gap equation 

2 



= 0, 



(25) 
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To determine the second-order phase transition between 
the normal and the superconducting phases, we take the 
limit Ai ^ in the gap equation where the supercon- 
ducting gap vanishes smoothly. The resulting equation 



ffi y^J tanh[/3c(|ek| 



m)/2] 
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IfkllMl 
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tanh[/3e(|£k| +m)/2] 
kkl +M 



kk||/x| 
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determines the critical temperature /3c — l/ksTc for a 
given coupling gi and chemical potential /i, see Fig. 2(a). 

We first note that for zero chemical potential, the 
Fermi level lies at the conical points of the Dirac cone. 
Due to the vanishing of the density of states, the system 
is quantum critical. This means that the system can un- 
dergo a phase transition even at zero temperature, where 
there are no thermal fluctuations. In the present case, it 
is a second-order phase transition driven purely by quan- 
tum fluctuations. The quantum critical point, which is 
found to be gi,c/J — 6.2 in this channel, separates the 
normal and the superconducting phases. On the other 
hand, for finite chemical potential the system ceases to 
be quantum critical. The usual BCS picture of Fermi sur- 
face instability is then recovered where an infinitesimal 
attractive interaction favors superconductivity. 

Secondly, for small chemical potential /x/J <C 1, we 
observe that the Tc-curvc displays a non-monotonous de- 
pendence on the coupling around the temperature re- 
gion kBTc/J ~ 0{1). This is due to the intrinsic na- 
ture of the inter-band pairing in this channel. Inter-band 
pairing generally requires more energy fluctuations, since 
there is an energy gap in the pairing between a state 
from the upper energy band and its time-reverse part- 
ner from the lower energy band. For vanishingly small 
chemical potential, wc arc left with two energy scales in 
the problem, namely the thermal energy and the cou- 
pling. Whenever thermal fluctuations are suppressed for 
UbT / J < 1, a stronger coupling is therefore required 
to promote the pairing. Thus, we see the increase in 
the coupling strength required to induce superconductiv- 
ity as thermal fluctuations are cut off aroimd the region 
fcsT/ J ^ 0(1). This feature is called a re-entrant behav- 
ior because even when the system is in the symmetry- 
broken (ordered) phase below the critical temperature, 
it can get back to the disordered phase by lowering fur- 
ther the temperature. This is only valid for a coupling 
around the critical value g\^c and a small chemical poten- 
tial. To conflrm the re-entrant behavior, we also investi- 
gate the free energy function at various temperatures for 
a flxed coupling and zero chemical potential, as shown in 
Fig. 2(b). The absolute minimum of the free energy in- 
deed returns to the origin where the system is disordered, 
as the temperature is lowered. 
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E. Extended Superconductivity 

Now, we want to study the extended superconducting 
channel. In this case, we set Ai = in the Hamiltonian 
(23). Furthermore, for the sake of simphcity, we assume 
the four-component RVB order to be real. We then arrive 
at the following properties for the superconducting order 
parameter: 



^2,k 



Im A2 k 



+ ^2(e;)A2(e™) cos(k • (e; - e„j)), 
E A2(e;)sin(k- e;). 



The quasi-particle spectra are then given by 



i?i,2,k = y|ek|2 + |A2,kP + /i2T2|ek|^(ImA2,k)2+A/2 
Upon cxtremizing the free energy function 

^ -F[A2(ei),A2(e2),A2(e3),A2(e4)] =0, (28) 



for / ~ 1,2,3,4, we obtain a set of four coupled gap 
equations 



Ir ^ 



/tanh(/3£;i,k/2) tanh (/3£;2,k/2) 



El, 



E- 



2,k 



A2(e;) + E A2(e,„) cos[k • (e/ - e,„)] 
tanh(/3£;i,k/2) ^ tanh (/3i;2,k/2) 

£^l,k £^2,k 

Ekl sin(k • e;) Im A2,k 



V(ImA2,k)2+M2 



for Z = 1, 2, 3, 4. We now take the system to be close to 
the phase transition, where the gap A2(e/) is small, and 
expand the right-hand side of the gap equations to lead- 
ing order in the gap to obtain the linearized gap equations 



m— 1 k 

tanh[/3e(|ek| - 



/ tanh[/3,(|£k| + |/i|)/2] 

I kkl + IH 



cos[k • e;] cos[k • e„ 



sinh(^cl^l) 



Ml cosh[;3e(|ek| + |mI)/2] cosh[/?e(|ek| - |m|)/2] 



X sin(k • e/) sin(k • e^) 



A2(em) 



(29) 



^ -1 

O 
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FIG. 3: Solutions to the gap equation with (a) d-wave symme- 
try and (b) p-wave symmetry at the phase transition for differ- 
ent chemical potentials: p/J — >■ (solid), /i/J = 1 (dashed), 
and /i/J = 2 (dotted). 



for I = 1,2,3,4. The linearized gap equations can also 
be written in the matrix form 



n = 



92 



fD B C B\ 
B D B C 
C B D B 

\b C B Dj 



where the four-vector f2 is defined by ff" = 
(A2(ei),A2(e2),A2(e3),A2(e4)), and D,B,C are ma- 
trix elements obtained from the coupled Eqs. (29). In 
fact, a numerical evaluation shows that the matrix ele- 
ment S ~ 0. It is clear now that the problem of de- 
termining the critical temperature Tc amounts to finding 
non-trivial solutions to the matrix equation. The eigen- 
values and their corresponding eigenvectors are classified 
as follows: 



C + D 
92 \-C + D 



(1,0, 1,0), (0,1, 0,1), d-wave 
(-1,0, 1,0), (0,-1, 0,1), p-wave 

(30) 
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with 



C + D 



1 V- 



2N 



tanh[/3e(|ek| + \n\)/2] 
kk| + |Ai| 



-C + D 



sinh(^c|/"l) 



2iV 



|/x|cosh[^e(|ek| + |Ml)/2] 
sm2(k+d) 



cosh[;3e(|ek|-|M|)/2]' 

where = {k^ ± ky)/2. The four eigenvectors come in 
two classes of symmetries, which reflect the underlying 
irreducible representation of the symmetry group of the 
square lattice. In Fig. 3, we solve for the critical temper- 
ature as a function of the coupling, for the two classes of 
eigenvectors. Again, a quantum critical behavior is ex- 
pected for zero chemical potential. However, in contrast 
to the local pairing in the previous section, this channel 
does not involve inter-band pairing and the system should 
not display re-entrant features. In Fig. 3(a) we show the 
Tc-curves for the rf-wave channel which behave as ex- 
pected. The quantum critical coupling is g2,c/J — 2.1, 
which is smaller than the critical local pairing coupling. 

On the other hand, the T^-curves for the p-wave chan- 
nel exhibit highly irregular features, see Fig 3(b). To 
understand that this is an artifact of the solution of the 
gap equation, let us analyze the free energy function in 
the p-wave channel more closely. By decreasing the tem- 
perature below the critical value given by the upper part 
of the Tc-curvc at a fixed coupling, the change in the free 
energy is shown in Fig 4(a). While the absolute minimum 
is shifted away from the origin, signaling that the system 
enters the symmetry-broken phase, a new local minimum 
is also being developed at the origin below the second 
critical temperature. Since the linearized gap equation 
is an expansion around the origin, the lower part of the 
Tc-curve seen in Fig 4(b) merely describes the develop- 
ment of the new local minimum. We thus conclude that 
there is no re-entrant behavior. However, the develop- 
ment of another local minimum at the origin gives rise to 
the possibility of a first-order phase transition where the 
Tc-curve develops a vertical slope. Looking at the free 
energy, we indeed find a first-order phase transition line 
starting at a tricritical point, where it meets the upper 
part of the second-order phase transition Tc-curve. 

Given the choices of pairing states with different sym- 
metries, we need to determine the most favorable one by 
comparing the free energy deep in the superconducting 
phase. Comparing the Tc-curves for the two channels in 
Fig. 3(a) and Fig. 3(b), we see that the rf-wave channel 
generally has a higher critical temperature for coupling 
92 /J 2.8. For 92/ J ?J 2.8, since the p-wave chan- 
nel can be favorable, we need to consider more general 
states with mixed symmetry. As already discussed in 
Sect. IIIB, the lack of parity symmetry for the pairing 
state in this channel allows for a coherent mixing of order 
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FIG. 4: (a) Evolution of the free cnergjr for the p- 
wave superconductivity at different temperatures ksT/J = 
3 (dotted), 1.5 (dashed), 0.1 (solid), at a fixed coupling 
92/ J = 5 and chemical potential n/J = 0. (b) Comparing 
the free energy for the p-wave channel (- A2/ ^2, 0, A2/ ^2, 0) 
(dotted), the d-wavc channel (A2/\/2, 0, A2/\/2, 0) (dashed), 
and the mixed channel (A2, 0,0,0) (solid) which yields the 
lowest free energy, for /it = 0, fcsT/J = 0.5, 92/ J = 



parameters with different symmetries. By taking the four 
eigenvectors as a basis which spans the space of the order 
parameter, we look for the vector which yields the low- 
est free energy. As shown in Fig. 4(b), wc find the state 
(1, 0, 0, 0) to be the most favorable (lowest free energy) in 
the superconducting phase and the corresponding order 
parameter reads 



A2,k = A2 



cos 



I sm 



(31) 



A comparison of the real-space configurations of the lo- 
cal and non-local pairing, see Figs. 7(a) and 7(b), shows 
that the latter leads to unconventional superconductiv- 
ity, where both the gauge and the €4,^ crystal lattice 
symmetries are spontaneously broken. 



F. Phase Diagram at Quantum Criticality 

Until now we treated the two superconducting channels 
separately. Here we consider the competition between the 
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FIG. 5: (color online) (a) Real-space configuration of the local 
pairing s-wave superconducting phase, (b) Real-space config- 
uration of the nearest-neighbor spin-singlet bonding with one 
no n- vanishing component. 



different types of superconductivity at the quantum crit- 
ical regime with zero chemical potential and a staggered 
flux (j) = TT. The discussion for general values of chemical 
potentials and flux values will be presented in the next 
section. 

We have shown that the most favorable supercon- 
ducting state in the separate channels are the constant 
gap Ai and the unconventional superconducting state 
^2(6;) — (1,0,0,0), respectively. We choose the latter, 
even though we have already noted that the extended 
superconducting channel has a much richer behavior 
around the region g2/J — 2.8. By taking into account 
the possibility of co-existence between the two supercon- 
ducting phases, we numerically minimize the free energy 
F(Ai, A2,k) containing both superconducting orders. By 
locating the global minimum in the free energy, we find an 
unpaired phase (normal phase) , an s-wave superconduct- 
ing phase and a non-local pairing superconducting phase 
in the phase diagram, see Fig. 6. The two superconduct- 
ing phases are separated by a first-order phase transition 
(dashed line) with no region of coexistence. A multi- 
critical point is identified at {gi/ J, 92/ J) = (6.2,2.1) at 
zero temperature, where the first-order phase transition 
line and the two second-order lines meet. 

It is important now to justify the validity of the mean- 
field approach to the superconductivity problem. The 
general criterion for the validity of the mean-field theory 
in the weak coupling regime is that the coupling strength 
must be smaller than the energy bandwidth 5E. For the 
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Unconventional 
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FIG. 6: (color online) Zero-temperature mean-field phase di- 
agram for the Hamiltonian (19) with zero chemical potential 
and a staggered flux tf) — tv. The two superconducting phases 
are separated by a flrst order phase transition line (dash line) 
between them, and a second order phase transition (full line) 
with the normal phase. 



system we study, the energy bandwidth is given by 



6E 




(32) 



The weak coupling regime is then given 171,2 < SE. Re- 
ferring to the phase diagram in Fig. 6 and a bandwidth 
SE ~ 5. 7 J, we find that both the 171 and 172 channels 
are approximately within the weak coupling regime, even 
though the 171 channel is closer to an intermediate cou- 
pling regime. Nonetheless, for coupling strength within 
the weak coupling regime, the three distinct many-body 
phases are already accessible. 



G. Phase Diagram Away from Quantum Criticality 

In this section, we consider arbitrary values of the 
chemical potential. In this case, rather than the chemical 
potential, the physical quantity that can be controlled in 
experiments is the particle density, or the fermion filling 
fraction (n), the average number of spin- 1/2 fermions per 
lattice site. Furthermore, comparing the two sets of solu- 
tion for the gap equations in the different channels at zero 
chemical potential, see Fig. 2(a) and Fig. 3(b), we note 
that the unconventional superconducting channel yields 
the highest transition temperature Tc- Thus, we shall 
consider only this channel with a varying fermion filling 
fraction. At the phase transition T — Tc, where the su- 
perconducting gap vanishes, the fermion filling fraction 
(n) can be determined quite easily by the non-interacting 
limit (ignoring the Hartree energy) 



\{n) 



II 



N ^ 



2N 



tanh 



IMI 



2fcBT, 



- tanh 



|ek| 



(33) 
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FIG. 7: (color online) Finite-temperature phase diagram in 
the unconventional superconducting channel for different dop- 
ing 5, for gi/J < 6.2. The shaded plane at a fixed coupling 
Qi/J ~ 1 shows a dome-like shape superconducting phase 
embedded in the normal phase. The latter interpolates Dirac 
liquid behavior at low-doping to a Fermi liquid behavior at 
high-doping. 

where N is the total number of sites. The quantity 6 
is conventionally called hole (particle) doping in solid 
state materials, since it measures the departure of the 
electronic density from the half-filling (particle-hole sym- 
metric) limit. By solving the doping equation (33) self- 
consistently with the linearized gap equation (29), we 
obtain the Tc-curve as a function of doping summarized 
in Fig. 7. 

We see that as soon as the system is tuned away from 
unit filling fraction {fj, — 0, S = 0), the system ceases 
to be quantum critical (The exponential tail in the Tc- 
curves extending to down to g2/J = is however not 
visible in the temperature scale of Fig 7). Secondly, the 
shift of the Tc-curves is not a monotonous function of 
the filling fraction. The Tc value increases as the filling 
fraction decreases from unity, but below a filhng fraction 
of approximately 0.3, the value decreases again (see 
the Tc-curves for different filling fractions in Fig. 7). 

To understand better the non-monotonous shift of the 
Tc-curves, it is illustrative to observe the evolution of the 
non-interacting Fermi surface. As shown in Fig. 8, as 
we tune the doping from zero to unity, the nature of the 
fermionic carriers changes from hole-like to particle-like 
at the doping 6 — 0.5, or the quarter-filling fraction. For 
S < 0.5, the Fermi wave vector kp increases in the two 
inequivalent Fermi pockets as the doping increases. For 
S > 0.5 the Fermi wave vector kp decreases instead as 
the doping continues to increase. Qualitatively speak- 
ing, the non-monotonous shift of the Tc-curve is based 
on the behavior of the Fermi wave vector. The BCS 
picture relates the transition temperature to the Fermi 
wave vector via an exponential function, ~ e""'^/'^^'"' , 
where a is proportional to the scattering length in an 
atomic system. This means that for an increasing Fermi 
wave vector in the doping region S G [0, 1/2], there will 
be an increasingly higher transition temperature, while 
for the decreasing Fermi wave vector in the doping region 



6 G [0, 1/2], the transition temperature decreases. 

Finally, by plotting temperature versus doping for a 
fixed value of the coupling strength gi, as shown in the 
shaded plane in Fig. 7, we find a dome-shaped unconven- 
tional superconducting phase at intermediate filling frac- 
tions, surrounded by the normal phase for fillings close 
to zero or unity, which we termed Dirac-liquid (Fermi- 
liquid) on the left (right) side of the phase diagram, where 
linear (quadratic) single-particle dispersion prevails. The 
dome-structure of the unconventional superconducting 
phase is similar to the phase diagram for high- Tc cuprates 
and heavy fermions. 

H. Away from Isotropic Dirac Cones 

So far we have fixed the staggered flux value to be 
(j) = TT. Except for the special values of </> = 27rz/, i/ S Z 
we do not expect a qualitative change in the mean-field 
results. The effect of the different staggered flux values 
is a change in the anisotropy of the Dirac cones. For a 
flnite doping, it results in anisotropic Fermi surfaces, see 
Fig. 9. 

IV. EXPERIMENTAL CONSIDERATIONS 

The superconductivity considered in this article arises 
for temperatures on the order of ten percent of the Fermi- 
temperature. This temperature range can be accessed in 
state of the art Bose- Fermi mixtures subjected to an opti- 
cal lattice. A possible experimental scenario is to employ 
the widely used rubidium-potassium system composed of 
a balanced mixture of fermionic ^"K-atoms prepared in 
the \F = 9/2, mp = -7/2) and \F = 9/2, mp = -9/2) 
Zceman components of the F — 9/2 ground state hypcr- 
flne level and bosonic ^^Rb atoms in the |f — 1, mp = 1) 
ground state [41]. The parameter U^p/Ubb may be ad- 
justed via its dependence on the well depth, while an s- 
wave Feshbach resonance around 202 Gauss can be used 
to tune Upp independently. In order to experimentally 
discriminate the two SC phases discussed here, one could 
search for a signature of their distinct gap functions in 
their momentum spectra [43] . Correlation measurements 
similar to the one described in Ref. [44] should be another 
powerful method to obtain information on the nature of 
the pairing. 

V. DISCUSSIONS AND CONCLUSIONS 

We considered here an ultracold Bose-Fermi mixture 
in a 2D square optical lattice subjected to an effec- 
tive staggered magnetic fleld, which exhibits a Dirac-like 
single-particle spectrum. The system is studied at low- 
temperatures, such that the bosons condense and me- 
diate a longer-range (nearest neighbor) attractive inter- 
action between the fermions. At half-fllling, the Dirac 
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(a) 




FIG. 8: (color online) The intersection of the shaded plane with the lower branch single-particle spectrum represents the 
topology of the Fermi pockets. The filling fraction determines the height of the shaded plane. From figures (a)-(c) the filling 
fraction is decreasing. 





(a) 



(b) 



FIG. 9: (color online) (a) For flux 4> ^ it, the Fermi pockets 
generally take the shape of "banana" due to the anisotropic 
Dirac cones, (b) At lower filling fraction, the resulting sets of 
nesting vector are also different from the 7r-flux case due to 
Fermi "squares" of different sizes. 



fermions experience both, local (Hubbard) and nearest 
neighbor attractive interactions, which can be tuned in- 
dependently. The zero-temperature mean field phase dia- 
gram exhibits a competition between a local s-wave and 
a non-local unconventional superconducting phases. It 
is interesting to compare the superconductivity occur- 
ring here to that of graphene-like systems [45-47]. In 
the square lattice with a staggered flux, equivalent Dirac 
cones are related to each other by time-reversal. In con- 
trast, in the graphene lattice, the time-reversal operation 
maps one Dirac cone to the other inequivalent Dirac cone. 
Thus, when BCS Cooper pairs are formed in the square 
geometry, no different flavors of Dirac fermions are in- 
volved. This is not the case for the hexagonal symmetry. 

At finite temperatures, the appearance of unconven- 
tional superconductivity below a dome reveals an intrigu- 
ing link to strongly correlated electronic materials. In 
addition, the evolution of the normal phase (surround- 
ing the superconducting dome) from a Dirac-liquid to a 
Fermi-liquid upon increasing doping is another essential 
feature of high- Tc cuprates and heavy fermions, not eas- 
ily captured by usual theoretical descriptions of the sys- 
tem. Finally, the appearance of "banana-shaped" Fermi- 



pockets for flux values 7^ tt, reminiscent of the ones 
observed in high-T^ cuprates by means of ARPES ex- 
periments [48] , adds to the number of puzzling similari- 
ties. The local anisotropics in the Fermi surface in Fig. 9 
could eventually lead to striped ground states, like those 
observed in the cuprates [49]. The multitude of similari- 
ties between the system considered here and the high-Tc 
superconductors suggests that also for high-Tc-niaterials 
a generalized Hubbard model with complex rather than 
real hopping coefficients might be the appropriate de- 
scription. This speculation presumes that some physical 
mechanism could be breaking the time-reversal symme- 
try already in the pseudogap phase of cuprates, which 
is in fact supported by recent observations [50]. A re- 
maining open question is then to identify which particle 
should be playing the role of the bosons in the cuprates, 
to mediate a longer-range interaction, which ultimately 
leads to unconventional superconductivity. In any case, 
cold atomic systems are already proving to be a fascinat- 
ing playground to develop our understanding of high-Tc 
superconductors. 
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VI. APPENDIX 

Here we provide the details of the Bogoliubov approx- 
imation procedure for bosonic operators to arrive at the 
mean-field action Eq. (9). The procedure amounts to 
selecting out the condensation mode and making the re- 
placement /3ko ~^ V^o + Po in the Hamiltonian. The 
Bogoliubov approximation amounts to keeping only the 
terms with fluctuation operators up to second order. For 
convenience, we deflne the various coupling constants as 
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Oi = Ubb/4:N, 02 = Ubf/N. For the terms containing only bosonic operators we have 



Hb = ^(i;k-M)4/3k + ai ^ M(ki,...,k4)4^4^^k3/9k4^(ki+k2-k3-k4) 

k ki,k2,k3,k4 



{Eo -n)+ 4ai7Vo 



7Vo(/3o + 4) + 2aiA^o' + I^(^^k - M + 8aiNo)plpu 



+aiiVo ^ M(ko, ko, k, -k)/3k/3-k + aiiVo ^ M(k, -k, ko, ko)/3^/3lk, 

k k 

where we have used M(ki,k2,ki,k2) = 2. And for the boson- fermion interaction term on the ^-sublattice we get 
Hbf,a = a2EE"r,.nf = a2 5]^Lr,.5]5k^4^e-*''---^5k./3k.e*-4 

re^ cr re^ o- I ki k2 ) 

re^ ^ ^ L k 

+ \/^Soe"^° '-E 3^^/316-^"--+ ^ <?k,5k./3l,/3k.e-'''--+*--j| 

k ki,k2 -I 

~ azNoN J2^'^ + 0'2^/NoN ^ n„{po + 4) + aiNh^ ^ ^^k + a2iVo Yl Yl """-'^ 

(T k reA f 



ik-r 



where we have introduced a mean-field density ha for the spin a fermions and thus, fluctuation terms of higher order, 
i.e., O (?^r,a•/3k/3k)^ Can be neglected. Performing the same procedure for the boson-fermion interaction term on the 
S-sublattice yields a similar expression Hbf,b- CoUecting aU the terms we get 

Hb + Hbf, a + Hbf,b 

= (eo-ij. + UBBno/2 + UbfJ2^A^o + \^(eo - m + UbbUq + UbfJ2^'^)(I^o + 4) 
+ X] (^k - M + 2f/BBno -h Ubf riajplPi, + ^UbbUq ^ ^(^o, ko, k, -k)/3k^-k 

k ^ <T ^ k 

+ ^UBBno Y ^i^' ^0, ko)/3k/3lk + Ubf no ^ X(nr,<T + rir+ei.a) 

k r£^ CT 

+f^^^\/^E|E E-r,.5S5ke-^(>^->^)- + X Xnr+e.,.e-^(''-'')-(-+-)|/3k 

+UBFMj:{j:j:^r,.909le^^^"-^y'^ + E E-r+e^.e^^-^-'^Hr+eoUt. 



(34) 



re^ cr 

r 
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